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Abstract 

Shrinkage algorithms are of great importance in almost every area of statistics due to the increasing impact 
of big data. Especially time series analysis benefits from efficient and rapid estimation techniques such as the 
lasso. However, currently lasso type estimators for autoregressive time series models still focus on models with 
homoscedastic residuals. Therefore, an iteratively reweighted adaptive lasso algorithm for the estimation of time 
series models under conditional heteroscedasticity is presented in a high-dimensional setting. The asymptotic 
behaviour of the resulting estimator is analysed. It is found that the proposed estimation procedure performs 
substantially better than its homoscedastic counterpart. A special case of the algorithm is suitable to compute 
the estimated multivariate AR-ARCH type models efficiently. Extensions to the model like periodic AR-ARCH, 
threshold AR-ARCH or ARMA-GARCH are discussed. Finally, different simulation results and applications to 
electricity market data and returns of metal prices are shown. 

Keywords: High-dimensional time series. Lasso, Autoregressive process. Conditional heteroscedasticity. 
Volatility, AR-ARCH 


1. Introduction 

High-dimensional shrinkage and parameter selection techniques are of increasing importance in statistics 
in the past years. In recent years, high-dimensional shrinkage and parameter selection techniques have been 
of increasing importance. In many statistical areas, lasso (least absolute shrinkage and selection operator) 
estimation methods, as introduced by Tibshirani (1996), are very popular. In time series analysis the influence 
of lasso type estimators is growing, especially as the asymptotic properties of stationary time series are usually 
very similar for stationary time series to the standard regression case, see e.g. Wang et al. (2007b), Nardi and 
Rinaldo (2011) and Yoon et al. (2013). Hence, given the lasso’s shrinkage properties, it is attractive for subset 
selection in autoregressive models. In big data settings, it provides an efficient estimation technique, see Hsu 
et al. (2008), Ren and Zhang (2010), and Ren et al. (2013) for more details. 

Unfortunately, almost the entire literature about £i-penalised least square estimation, like the lasso, deals 
with homoscedastic models. The case of heteroscedasticity and conditional heteroscedasticity is rarely has 
rarely been covered so far. Recently, Medeiros and Mendes (2012) showed that the adaptive lasso estimator is 
consistent and asymptotically normal under very weak assumptions. They proved that the consistency and the 
asymptotic normality hold if the residuals are described by a weak white noise process. This includes the case 
of conditional heteroscedastic ARCH and GARCH-type residuals. Nevertheless, their classical lasso approach 
does not make use of the structure of the conditional heteroscedasticity within the residuals. Without going into 
detail, it is clear that the estimators might be improved if the structure of the conditional heteroscedasticity in 
the data is used. Furthermore, Yoon et al. (2013) analysed the lasso estimator in an autoregressive regression 
model. Additionally, they formulated the lasso problem in a time series setting with ARCH errors. However, 
they did not provide a solution to the estimation problem and left this for future research. 

Recently, Wagener and Dette (2012) and Wagener and Dette (2013) analysed the properties of weighted 
lasso-type estimators in a classical heteroscedastic regression setting. They showed that their estimators are 
consistent and asymptotically normal. In addition, their estimators perform significantly better than their 
homoscedastic counterpart. Their results, conditioned on the covariates, can be used to construct a reweighted 
estimator that also works in time series settings. 

We derive an iteratively reweighted adaptive lasso algorithm that addresses the above mentioned problems. 
It enables the estimation of high-dimensional sparse time series models under conditional heteroscedasticity. 
We assume a regression structure which is satisfied by the majority of the important time series processes and 
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which admits fast estimation methods. The computational complexity of the algorithm is essentially the same 
as the coordinate descent algorithm of Friedman et al. (2007). This very fast estimation method for convex 
penalised models, such as the given £i situation, can be applied to the iteratively reweighted adaptive lasso 
algorithm. 

The algorithm is based on the results of Wagener and Dette (2013), as their results can be generalised to 
models with conditional heteroscedasticity. The sign consistency and asymptotic normality for the proposed 
estimator is adduced. Furthermore, a general high-dimensional setting, where in which the underlying process 
might have an infinite amount of parameters, is considered. Note that all the time series results hold in a 
classical regression setting as well. 

However, we restrict ourself to £i-penalised regressions as they are popular in time series settings (see e.g. 
Wang et al. (2007b), Nardi and Rinaldo (2011) and Yoon et al. (2013)). In general, other .^g-penalty could also 
be considered, e.g. the penalty. The £2 penalty, which gives the ridge regression, is suitable for shrinkage 
as well, but does not allow for sparsity. However in £q-penalised regression, the case g = 1 is the greatest case 
of practical intereststill allowing for sparsity. This sparsity property can be used in applications to select the 
required tuning parameter based on information criteria that are popular in time series analysis. 

The general problem ist stated in section 2. In section 3, we motivate and provide the estimation algorithm. 
Subsequently, the asymptotics are discussed in in section 4.In Section 5, an application to multivariate AR- 
ARCH type processes is considered. This includes several extensions such as periodic AR-ARCH, AR-ARCH 
with structural breaks, threshold AR-ARCH and ARMA-GARCH models. The section 6 shows simulation 
which underline the results given above. It provides evidence that incorporating the heteroscedasticity in a 
high-dimensional setting is more important than in low dimensional problems. Finally, we consider the proposed 
algorithm as a model for the electricity market and metal prices returns data. A two-dimensional AR-ARCH 
type model is used in both applications to the hourly data, in the first one to electricity price and load data 
and in the second one to gold and silver price returns. 

2. The considered time series model 

The considered model is basically similar to the one used by Yoon et al. (2013) or Medeiros and Mendes 
(2012). Let {Yt)t^z be the considered causal univariate time series. We assume that it follows the linear equation 

Yt = X^,t(3l+et, (1) 

where Xao,t = (ATi,*, Ar 2 ,t,...) is a possibly infinite vector of covariates of weakly stationary processes {Xi^t)t&, 
(et)tgz is an error process, and the parameter vector is = (/3°,/32, ■ • 0^ with l/^?! < o®- The covariates 

can also contain lagged versions of Yj, which allows flexible modelling of autoregressive processes. 

A simple example of a process that helps for understanding this paper is an invertable seasonal MA(1) 
process. In particular, the AR(oo) representation of a seasonal MA(1) with seasonality 2 is useful. It is given 
by Yt = Et - 9 et -2 = dYt -2 + d‘^Yt -4 + O^Yt-e + ... + Et, choosing Xoo.t = (Yt_i, Yt_ 2 ,...) with /3^ = 
(0,0,0,0^, 0,0^, 0,...)'. The error process {Et)t^i is assumed to follow a zero mean process with Ct being 
uncorrelated to the covariates X^o^t- Hence we require E(et) = 0 and Cov(et, Yi_t) = 0 for all f S N. Moreover, 
we assume that Et is a weak white noise process, such that 

Et = (JtZt where ct* = 5 ( 0 :^; Too,t) and {Zt)t& is i-i-d. with E{Zt) = 0 and Var(Zt) = I. (2) 

Here, g is a positive function, Xoo.t = (Ai,t, ^ 2 ,* j...) is a possibly infinite vector of covariates of weakly stationary 
processes {Li^t)tEi.i and cx.^ = (aj, a®,.. ■)' is a parameter vector. Similarly to the covariates Xao,t in (1), Too.t 
can also include lags of at or Et- This allows for a huge class of popular conditional variance models, like ARCH 
or GARCH type models. Choosing 

g{cx°^; Loo,t) = g((ao, ai, ■ ■ •); (et-i, crt_i, 0, ...)) = + aie?_i +a 20 't_i 

leads to the very popular GARCH(1,I) process. Note that the introduced setting is more general than the 
conditional heteroscedastic problem stated by Yoon et al. (2013), who mentioned only ARCH errors. 

For the following we assume that the time points 1 to n are observable for Yt. Thus, we denote by 

/ \ / Yi 1 • • • All \ ( \ 

Y ji I I , X I I , /3^ I I , and Yyj JYyj/3^ 

\ Yn j \ A^n,l * * ' Xji^p^ j y ^’Pn. j 

the response vector the nxpn matrix of the covariates X^ the parameter vector /3„ and the corresponding 
errors e„. Furthermore let Xi ,be the rows of X„. 
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Since we deal with a high-dimensional setting we are interested in situations where the number of possible 
parameters increases with sample size n. Therefore, denote /3° = (/3 °,..., Pp^)' the restriction of /3^ to its 
first pn coordinates. Due to \l^i \ < oo it follows for e° = = Yn — Xnf3n that there is 

a positive decreasing sequence (Cn)n with 0 such that lim„_>.oo T’(maxi<(<„ |e° ^ — £(| < C„) —1 holds. 

Thus, for a sufficiently large n we can approximate Yt by Xn,t0li arbitrarily well. 

However, under the assumption of sparsity, meaning that only some of the regressors attribute significantly to 
the model, we can conclude that only (/„ of the Pn parameters are non-zero. Hence, there are parameters 

that are exactly zero. Without loss of generality we assume that and ,3° are arranged so that the first 
components of /3° are non-zero, whereas the following are zero. Obviously we have /3° = (/3°,..., , 0,..., 0)' = 

(/3° (1)', o')'. This arrangement of the non-zero components is only used to simplify the notation, it is especially 
not required by the estimation procedure. Additionally we introduce the naive partitioning of X^ and /3^, in 
such a manner that /3„ = (/3„(1)',/3„(2)')', = (X„(l), X„(2)) and = (X„,t(l)', A:„,t(2)')' holds. 

Subsequently, we focus on the estimation of /3°, for which we will utilize a lasso-based approach for 13^. 
Henceforth, we achieve never a direct estimate for /3)^, but we can approximate it by (/3)j, O')'. 

3. Estimation algorithm 

The proposed algorithm is based on the classical iteratively reweighted least squares procedure. An example 
for an application of it to time series analysis can be found e.g. in Mak et al. (1997). However, similar approaches 
are not popular in time series modelling, as there are usually better alternatives if the number of parameters 
is small. In that case, we can simply perform an estimation of the joint likelihood function of (1), see e.g. 
Bardet et al. (2009). But when facing a high-dimensional problem, it is almost impossible to maximise the non¬ 
linear loss function with many parameters. In contrast, our algorithm can be based on the coordinate descent 
lasso estimation technique as suggested by Friedman et al. (2007) which provides a feasible and fast estimation 
technique. Other techniques, like the LARS algorithm introduced by Efron et al. (2004) which provides the full 
lasso solution path, can be used as well. 

For motivating the proposed algorithm, we divide equation (1) by its volatility, resp. conditional standard 
deviation at- Thus, we obtain 

= + (3) 

where Yj = and Xao,t = ~A^oo,t- Here, the noise Zt is homoscedastic with variance 1. Hence, if 
the volatility at of the process Yt is known, we can simply apply common lasso time series techniques under 
homoscedasticity. Unfortunately, this is never the case in practice. The basic idea is now to replace at by a 
suitable estimator at , which allows us to perform a lasso estimate on a homoscedastic time series as in equation 
(3). 

For estimating ARMA-GARCH processes, practitioners sometimes use a multi-step estimator. This estima¬ 
tion technique involves computing ARMA parameters in a homoscedastic setting first and then use the resulting 
estimated residuals are used to estimate the GARCH part in a second step, see e.g. Mak et al. (1997) or Ling 
(2007). We will apply a similar step-wise estimation technique here. 

In general, we have no a priori information about at, hence we should assume homoscedasticity in a first 
estimation step. We start with the estimation of the regression parameters /3°, resp. /3^, and obtain the 
residuals £„j,..., £„_„. We use the residuals to estimate the conditional variance parameters and thus 
(cti, ... ,CT„) by (fTn.i, ■ • ■ ,^n,n) afterwards. Afterwards, we reweight model (1) by a^^ to get a homoscedastic 
model version which we utilise in order to reestimate ,9° again. We can use this new estimate of /3° to repeat 
this procedure. Thus, we will end up in an iterative algorithm that hopefully converges in some sense to /3°, 
resp. with increasing sample size n. 

We use an adaptive weighted lasso estimator to estimate /3° within each iteration step. It is given by 

n / Pn X 2 Pn 

/3n.lasso(^n,'Mn,'U;n) = arg min ^ J Yj -^Xt,iPij + Xn^Vn,j\l3j\ 

^ t=l ^ i=l ' j=l 


or in vector notation 


(^nMsso{K,Vn,Wn) = argmin(- X„/3)'W^(Y„ - X„/3) + 

P 

where = diag(i(;„), Wn = (wn,i, ■ ■ ■, Wn,n) are the heteroscedasticity weights, Vn = {vn,i, ■ ■ ■, Vn,p„) are the 
penalty weights and A„ is a penalty tuning parameter. As described above, in the iteratively reweighted adaptive 
lasso algorithm we have the special choice Wn = {wn,i, ■ ■ ■ ,Wn,n) = for the heteroscedasticity 

weights within each iteration step. We require Wn = 1 for the homoscedatic initial step. 

Like Zou (2006) we consider, for the tuning parameter the choice = /3n Init fo^' some r > 0 and some 
initial parameter estimate f3^ With t = 0 we obtain = 1 which is the usual lasso estimator. Obviously, 
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there is no initial estimator required in this case. However, we consider the case of r = 0 and the adaptive lasso 
approach for our practical application, as they resulted in different perfomances. 

The selection of the tuning parameters A„ and r such as the choice of the initial estimate /3„ is crucial 
for the application and might demand some computational cost. We discuss this issue in more detail at the end 
of the next section. 

Subsequently, we denote Yn) as a known plug-in estimator for , which is the projection 

of to its first coordinates. We denote as restriction of g that corresponds to q:°. Thus, gn is defined 
such that is restricted to and is a restriction of = {Lao,t)t& to its first mn{ln) coordinates. 
Similarly, let X„, y„) be an estimator for {L^ i, ■ ■ ■, L^ ^Y- 

For example, if St follows a GARCH(1,1) process we receive at = gn{anT ^n,t) for all n S N, where q:° = 
(ao,ai,a 2 ) with = 3 and L^ t — (st-i, o't-i) with mn{ln) = 2 for all n € N. This is similarly feasible 
for every variance model with a finite amount of parameters. However, if at follows an infinite parameterised 
process, e.g. through an ARCH(oo) process. In and m„(Z„) should tend to infinity as n —>■ oo. 

The estimation scheme of the described iteratively reweighted adaptive lasso algorithm is given by: 

1. initialise A„ > 0, ti„(r) = (n„,i(r),... ,u„,p„(r)) = with r > 0 and = 1 , A: = 1 

2. estimate by weighted lasso: /3[f' = = /3„_iasso(A„, u„(r), 

3. estimate the conditional variance model: ctn^ = X„, y„) and y„) 

4. compute new weights ..., Wn}n) with lajf ^ = gn{an \ L^nl)~^ 

5. if the stopping criterion is not met, A; = fc -|- 1 and back to 2. otherwise, return estimate and 
volatilities = 9 n{ctn \ L^n^t) 

We can summarise that we have to specify the tuning parameter A„, the initial estimator f3n init with an 
inital value of r, and the initial heteroscedasticity weights To reduce the computation time it can be 

convenient in practice to choose r = 0 (lasso) or r = 1 (almost non-negative garotte). 

The stopping criterion in step 5 has to be chosen as well, such that the algorithm eventually stops. A 

\k] [fel 

plausible stopping criterion should measure the convergence of Wn , resp. crk . We suggest to stop the algorithm 
if \\aYY — < e for a selected vector norm |j • |j and some small e > 0. Nevertheless, in our simulation 

study, we realised that the difference in the later steps are marginal, so that stopping at A: = 2 or A: = 3 seems 
to be reasonable for practice. This will be underlined by the asymptotics of the algorithm as analysed below; it 
can be shown that, under certain conditions. A; = 2 is sufficient to get an optimal estimator if n is large. 


4. Asymptotics of the algorithm 


For the general convergence analysis it is clear that the asymptotic of the estimator will strongly depend 


on the (cond.) heteroscedasticity models (2) (esp. the formula for g) such as on the linked estimators and 
Ln- Despite that strong dependence we are able to prove sign consistency as introduced by Zhao and Yu (2006) 
and asymptotic normality of the non-vanishing components of /3[f^ in a time series framework. 

If we assume that the number of parameters does not depend on the sample size n, then we could make 
use of the results from Wagener and Dette (2012) to obtain asymptotic properties, as they prove sign consistency 
and asymptotic normality under some conditions for the weighted adaptive lasso estimator. 

The case where the number of parameters increases with n is analysed euivalently in a regression frame¬ 
work by Wagener and Dette (2013), but only for the adaptive lasso case with r = 1. They basically achieve the 
same asymptotic behaviour as for the fixed case, but it is clear that the conditions are more complicated 
compared to those of Wagener and Dette (2012). 

In the following we will introduce several assumptions, which allow us to generalise the results of Wagener 
and Dette (2013). One crucial point is the assumption that the process Yt can be parameterised by infinitely 
many parameters, so that the error term = Yn — Xnfi^, based on the restriction /3” of the true parameter 
vector /3|^, is not identical to the true error restriction n- contrast to the term is in general 

correlated. This has to be taken into account for the proof concerning the asymptotic behaviour. 

~[k] r. ~ [fc] 

For the asymptotic properties we introduce a few more notations. Let and Y^ = 

where = diag(u;lf^). Let denote the true volatility matrix and its estimate 


in the A:-th iteration. Additionally, we introduce F 


[fe] —[fc] '—[fc] ~[i] 

= ^(-X'„ ) as the scaled Gramian, where r„ = r„ = 
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is the unsealed Gramian. Furthermore, let and r„ denote the weight matrix and the Gramian 

that correspond to the true matrix The submatrices to /3° (1) are denoted by r„ (1), r„(l), and rjj(l). 

Similarly to Wagener and Dette (2013), we require the following additional assumptions, which we extended 
to carry out our proof: 

(a) The process (Yj, ..., is weakly stationary with zero mean for all m € N. 

(b) The covariates are standardised so that = 1 for alH € Z and i G N. 

(c) For the sequence of covariates of a fixed t there is a positive sequence ('&n)ne'N such that 

max ||X„,t(l)||2 = C>p('d„v^)- 

l<t<n 


(d) For the minimum of the absolute non-zero parameters bn = min{|/3° (1)|} and the initial estimator (3n init 
there exists a constant 6 > 0 so that 

lim P(6min{|^„ i„it(l)r} < 6„) =0. 

n—>-oo ’ 

(e) There exists a positive sequence (r„)„gN with —>■ oo such that 

lim P(max{|/3„ i„it(2)|^} < = 0. 

n—>-oo ’ 

(f) There are constants 0 < Ao.min < Ao,max and 0 < Ai^min such that the eigenvalues satisfy 

P{\o ,min ^ '^min (r„(i)) < •^max n (1)) < Ao ,m.ax) ^ Ij 


p 


^Ai^min ^ Aniin(Tyj (1)) ^ Amax(r'^ 


1 , 


for n —> 00 . 

(g) There is a positive constant tTmin such that 

0 ^ f^min ^ (*0:71 (/^n; ; ^n) j .^n,£ (/^n ’ ^n)) 

for all n > with G N, t G {1,..., n} and /3„ in an open neighbourhood of ,9° . 

(h) The volatilities have afinite fourth moment, so E(tT^) = E(g(Q:|^, Loo,tY) < 00 for all t. 

(i) For all n G N the estimator and Ln are consistent for and L^^i ,..., additionally 

| 5 (a^, Xoo.t)-' - y„), in.£(/ 3 °; AC„, = Op{^) 


for some {hn)neK with /i„n 2 Q as n —>■ 00. 
(j) It holds for A„, d„, p„, g„, 6 „, r„, and /i„ that 

log(ra)^<'^^> log(gn)3 Q 


y/nb-n 






(iib ^ 

(iv) 


sAJlog(n)i{‘*=i>] 


^ 0 


(v) 


hny/n 

XtiTtl 


(") ^ 


(vii) 




(viii) 

' ' \/n 


0 
0 
0 
0 


as n —>■ 00 . 

(k) There are positive constants Ci, C 2 and d with 1 < d < 2 such that 

P(|et| > x) < Cl exp{-C 2 x‘’‘). 

(k’) It holds for A„, p„, (?„ and that q g^g ^ 
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Assumption (a) is standard in a time series setting, (b) is the scaling that is required in a lasso framework, 
(d) and (e) are usual assumptions in an adaptive lasso setting (see e.g. Zou (2006) or Huang et al. (2008)). 
(f) gives bounds for the weighted and unweighted Gramian. (g), (h) and (i) postulate properties required for 
the heteroscedasticity in the model, (j) states some convergence properties that make restrictions to the grow 
behaviour within the model, especially the number of parameters Pn and the number of relevant parameters g„. 
(k) makes a statement about the tails of the errors. 

Using the assumptions above we can prove sign consistency and asymptotic normality. 

Theorem 1. Under conditions (a) to (j), where either (k) or (k’) holds, it holds for all k > 1 that 

lim P (sign(/3W) = sign(/3o)) = 1. 

Moreover it holds for S with ||^n ||2 = 1 that 

Vns„(fc)-1^; (/3W(1) - /3n(1)) ^ N{0, 1) 

in distribution, where s^(l) = ^(j(r„(l))“^^„ and s‘^{k) = 'C(i(r„(l))“^^„ for k>2. 

The proof is given in the appendix. Note that the variance s'^{k) for /c > 2 is substantially smaller than 
s^(l). Hence the estimator has minimal asymptotic variance for all k >2. 

Due to the general formulation of the theorem assumption (j) contains several assumptions on problem 
characterizing sequences. The convergence rate of the volatility model is relevant as well. If we have 

that hn — Op{l) (e.g. the variance model is asymptotic normal) then the three conditions involving are 
automatically satisfied by the other conditions. This reduces the relevant conditions in (j) a lot. 

There is one condition in assumption (j) involving dn that is given through assumption (c). As it holds 
that maxi<t<„ ||X„_((1)||2 = Op{'dnytqf) it characterises the structure of regressors. Obviously it holds that 
'dn = tDp{l) if /3^ contains only a finite amount of non-zero parameters, so —>■ c for some c G N as 

n —)■ oo. However, there are many other situations where = Op{l) holds. For example, if we have that 
^oo,t(l) is stationary. In the example above, where Yt follows a seasonal moving average process the process, 
-^oo.t(l) = {Yt- 2 ,Yt- 4 ,Yt-e ,...) is stationary. 

Furthermore, there is the option of (k) or (k’) in the theorem, (k) restricts the residuals to have an exponential 
decay in the tail, like the normal or the Laplace distribution. However, this can be replaced by the stronger 
condition (k’) in the theorem. In this situation, polynomially decaying tails in the residuals are possible. Here 
a specification of the constant d in (j) is not required, as (k’) implies directly (j) (iv), which means that (j) (i) 
is not used in this case. More details are given in the proof. 

As discussed in Wagener and Dette (2013) the assumption (k) or (k’) has an impact on the maximal possible 
growth of the amount of parameters in the estimation. There are situations where under assumption (k) 
can grow with every polynomial order, even slow exponential growth is possible. In contrast, given assumption 
(k’) this is impossible. Here Wagener and Dette (2013) argued that sign consistency is possible for rates that 
increase slightly faster than linearly, such as ~ n log(n), but not for polynomial rates like Pn for some 

i5 > 0. Wagener and Dette (2013) do not discuss this case for the asymptotic normality. In this situation, we 
can get an optimal rate of for the number of relevant parameters (having ~ 1, ^ ^ 1 and 

dn ~ 1), when we have a polynomial growth for p„. 

The quite general formulation in (i) can be replaced by a more precise assumption when a variance model 
is specified. For example, if we have a finite dimensional conditional variance model where q:„ is asymptotic 
normal, i.e. converges with rate of n~"i, and (3„ i— gn{cin{f3^, Xn,Yn), Ln^tif^n'y ^n,Yn)) is twice continously 
differentiable with uniformly bounded derivatives, then (i) can be satisfied by = c under some regularity 
conditions on and or its estimated counterparts S„ and If in contrast is increasing we will usually 
tend to get worse rates for /i„. 

In empirical applications, practitioners often just want to apply a lasso type algorithm without caring much 
about the chosen size of n and p„. They tend to stick all available n and Pn into their model as long as it is 
computational feasible. However, usually it is feasible to validate the convergence assumptions in (j) at least 
partially. Therefore, we have to estimate the model for several sample sizes n and a specified growth rate 
for Pn and A„. As we can observe the estimated values for of the model we can get clear indications for 
the asymptotic convergence properties. This also helps to find the optimal tuning parameter A„. The tail 
assumption (k) can be checked using log-density plots and related tests. The moment restriction (h) to the 
volatilities can be validated using tail-index estimation techniques, like the Hill estimator. 

Note that in the algorithm A„^is assumed to be the same in every iteration. It is clear that if we have two 
different sequences (A„)„gN and (A„)„gN that satisfy the assumptions of the theorem, we can use them both in 
the algorithm. For example we can use (A„)„gN for the first iteration and (A„)„gN for the subsequent iterations. 
This might help in practice to achieve better finite sample results. 
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For finding the optimal tuning parameters we suggest to use common time series methods that are based 
on information criteria. Zou et al. (2007), Wang et al. (2007b), Zhang et al. (2010) and Nardi and Rinaldo 
(2011) analyse information criteria in the lasso and adaptive lasso time series framework. Possible options for 
this information criteria are the Akaike information criterion (AIC), Bayes information criterion (BIG) or a 
cross-validation based criterion. Here, it is worth mentioning that Kim et al. (2012) discusses the generalised 
information criterion (GIG) in a classical homoscedastic lasso framework where the amount of parameters 
depends on n. They establish that under some regularity conditions the GIG can be chosen so that a consistent 
model selection is possible. 

For the initial estimate /3„ that is required for the penalty weights there are different options available. 
The simplest is the OLS estimator, which is available if < n. Another alternatives are the lasso (r = 0), 
elastic net or ridge regression estimator, see e.g. Zou and Hastie (2005). Remember that we require an initial 
estimate I3n,init the adaptive lasso case if r > 0. 

Note that Wagener and Dette (2013) described a setting with two initial estimators. One for the adaptive 
lasso weights as we do, and another one for the weight matrix VF„. The first estimator corresponds to our 
jjjjj, whereas the second inital estimator is not required, as we can initialise the volatility weight matrix Wn 
by the homoscedastic setting. A similar result was achieved by Wagener and Dette (2013) who showed that the 
homoscedastic estimator can be used as initial estimator in their setting. 

5. Applications to AR-ARCH type models 

In the introduction we mentioned that one of the largest fields of application might be the estimation of 
high-dimensional AR-ARGH type processes. Therefore, we discuss a standard multivariate AR-ARGH model 
in detail. Afterwards, we briefly deal with several extensions, the periodic AR-ARGH model, change point 
AR-ARGH models, threshold AR-ARGH models, interaction models and ARMA-GARGH models. 

Let Yt = (Fi,t,..., Yd,ty be a d-dimensional multivariate process and I? = {1,..., d}. 

5.1. AR-ARCH model 

The multivariate AR model is given by 

4>i,j,kYj^t-k + £i,t (4) 

jev keii.j 

for i S D, where 4>ij^k are non-zero autoregressive coefficients, A j- are the index sets of the corresponding relevant 
lags and Si^t is the error term. The error processes (£i,t)tGZ follow the same conditional variance structure as in 
(2), so where and (^j,t)tGZ is i-i.d. with E(Zi_t) = 0 and Var(Zi_t) = I. 

Now, we define the representation (4) that matches the general representation (1) by 

Yi^t = + £i,t 

for i gH where the parameter vector /3j = (<('i,i,fc)fcG/i,ij • ■ ■ j i4>i,d,k)keii d) ^ad the corresponding regressor 
matrix Xi^t = (1, {Xi^i^t-k)keii,ij • ■ ■ j {Yi^d,t-k)keii,d)- Note that this definition of is only well defined if all 
Ii,j for j G TA are finite, if one index set is infinite we have to consider another enumeration, but everything 
holds in the same way. 

Furthermore, we assume that e* = ■ ■ ,£d,t)' follows an ARGH type model. In detail we consider a 

multivariate power-ARGH process which generalises the common multivariate ARGH process slightly. Recently, 
Francq and Zakoi'an (2013) discussed the estimation of such power-ARGH(oo) processes and showed applications 
to finance. It is given by 

“ OtijO T 'y 'y Ol-i,j,k\£j,t—k\ b (5) 

keJij 

with Jij as index set and 6i as power of the corresponding at. The parameters satisfy the positivity restriction, 
so atfi > 0 and atj^k > 0- Moreover we require that the Si’s absolute moment exists. Obviously, we 

have 

/ \ l/'Si 

— f 4- y ( y ( —fc I j 

j&Vk&Jij ^ 

where oli = (a^^o, (ai,i,fc)/ce, • ■ •, (ai.<i.fe)fcGj^.d) and Li = ((£i^t_fe)fcgj, ^,..., (ed,t_fe)fcgj, ^). Similarly as for 
f3i^ oti is only well defined if all Jij for j G T> are finite. Otherwise we have to consider another enumeration. The 
case Si = 2 leads to the well known ARGH process which turns into a multivariate ARGH(p) if Ji j = {1,... ,p}. 

For estimating the ARGH part parameters we will make use of a recursion that holds for the residuals. This 
is given by 

1 ^ 2 ,— cKjQ 4- y ( y ( fcI 4- Ui^t 

j&V kGJi.j 


(6) 
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where 5^,0 = liOtifi, = hotij^k and Ui^t = - li) with 7 ^ = ^i{5i) = E\Zi^t\^'. Here, Ui^t is a weak 

white noise process with = 0. The fitted values ctJ of equation ( 6 ) are proportional to the crj up to the 

constant 7 ^. As 7 ^ is the absolute moment of Z^ t, it holds that 7 ^ = 2, if = 2. If (5^ = 1 and follows 

a normal distribution 7 ^ it is V2tt~'^ « 0.798. If £74 exhibits e.g. a standardised t-distribution we will observe 
larger first absolute moments 74 . 

Clearly, the true index sets Aj and Jij are unknown in practice. Thus we fix some index sets Iij{n) and 
for the estimation that can depend on the underlying sample size n. If the true index sets Aj and Jij 
are finite, then the choices Iij{n) = {!,..., max(A,j)} and Ji^j(n) = {!,..., max(Aj)} are obvious. If A,j and 
Ji^j are infinite, Ii^j{n) and jTij{n) should be chosen so that they are monotonically increasing in the sense that 
Iij{n - I) C I*j (n) and Jij{n - I) C Jij{n) with UnGN^hi(”’) = ^ ^.nd UnGN = N- The size of Iij{n) 

and Jij{n) is directly related to the size of the estimated parameters pi^n for „ and for It holds 

that Pi n = 1 + = 1 + Horc, 13i n and oti^n are the restrictions of (3i and oci 

to their first and li^n coordinates. 

For the estimation of /3j and „ we can apply the iteratively reweighted adaptive lasso algorithm as 
described in the previous section. However, we have to specify an estimation method for the variance part. In 
particular we require the estimators and Li, or more precisely their restrictions and Li^n to its and 
krii^niU^n) coordinates. For Li^n^f^i n^ i,n) we have the estimator 

Li^n^t — Li^n,t{^(3 in: ^ i,t) — \^i,t ^'i,n,t(3in\ 

which provides an estimate for and For the estimation of we suggest to minimise the 

problem 

A 4 4Q;j|j2, (7) 

where Ai^t = (1, iLi^n,t-k)keJi.i, ■ ■ ■, iLd,n,t-k)keJi,d): which corresponds to the plug-in version of equation ( 6 ). 
For the estimation of (7) a common non-negative least squares (NNLS) estimation technique can be considered. 
If the variance equation is high-dimensional approaches like the positive lasso are suitable as well. Hence high¬ 
dimensional lasso type algorithms with positivity constraint can be applied for the parameter estimation. But as 
the residuals in ( 6 ) only follow a weak white noise process, there are more advanced results for the asymptotic of 
this procedure required For the non-restricted adaptive lasso Medeiros and Mendes (2012) show sign consistency 
and asymptotic normality under certain conditions for such a situation with a weakly stationary error process. 

However, the simple NNLS estimation procedure can act as a shrinkage procedure as well, as some parameters 
can be estimated to be 0. This well known sparsity effect of NNLS settings was recently analysed by Meinshausen 
et al. (2013) and Slawski et al. (2013). Slawski et al. (2013) provided evidence that the NNLS approach is 
potentially superior to the positive lasso. We use the NNLS algorithm for the computational applications as 
described by Lawson and Hanson (1995). 

5.2. Periodic AR-ARCH model 

Another class of models where we can apply the proposed estimation technique is the class of periodic AR- 
ARCH models. Here, we assume a model as described above, but all parameters are allowed to vary periodically 
over time. This is very suitable for modelling seasonal effects in high-dimensional data. 

Thus, the model for the conditional mean equation is given by 

u,=*,„(()+i: i: 4^i,j,k{l')Zj t—k T ( 8 ) 

jev keii.j 

and for the conditional variance equation 

A ,4 — 3- ^ ^ ^ Cli,j,k{t)\Sjd-k\ '■ (9) 

keJij 

As mentioned, the time dependent parameters vary periodically over time. Assuming a periodicity of S we 
have (j^i ^it) = 0 ,/(f 5 ~ ^ 41 ^ 2 , 0 ,/(f 1 and c^i j }^{t) = 

^i, 3 ,k,i{t)o:i,j,k,i, where Bi^^i and Bij^k,i are 5'-periodic basis functions. 

Note that the processes is in general not weakly stationary anymore. However, they are periodically weakly 
stationary (also known as weakly cyclostationary). So if S' e N then the subsequences {Yst-\-s)tGZ follow a 
weakly stationary process. For more details see e.g. Aknouche and Al-Eid (2012). 

As choice for the periodic basis functions, periodic indicator functions are suitable if S is small, the parameter 
space will be blown up by a factor of S. If S is large, a Fourier approximation, periodic B-splines or periodic 
wavelets might be a good choice as basis to keep the parameter space reasonable. 

As mentioned, the process Y 4 is not stationary in general, so the asymptotic theory given above can not 
be applied. Nevertheless, a similar theorem is likely to hold true for periodic stationary processes. In order to 
proof this statement one would have to focus on the level of the mentioned weakly stationary subsequences, 
similarly as in Ziel (2015). The estimation procedure can be then performed as in the AR-ARCH model part. 
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5.3. AR-ARCH with structural breaks 

Another field of possible applications is the one of change point models, i.e. models where we have at least 
one structural break. Here, the basic model is a time-varying AR-ARCH model as defined in equations (8) 
and (9) for the periodic AR-ARCH model. The basis functions are defined so that they can capture structural 
breaks instead of periodic effects. The resulting model is of the same structure as the change point model used 
by Chan et al. (2013). If we have a priori information about the change point we can take this into account. 
If we have no information, some clever segmentation of the time should be considered. One option is to allow 
a change in every parameter (especially (jjifi) and at every time point. This can be handled by choosing n 
basis functions for each parameter so that they build a triangular matrix. The resulting model is a special 
case of the so called fused lasso (see e.g. Tibshirani et al. (2005)) and suitable for change point analysis. This 
particular mentioned approach of modelling change points is analysed in Levy-leduc and Harchaoui (2008) and 
Harchaoui and Levy-Leduc (2010). However, this increases the parameter space enormously, in every case we 
receive Pn > ri. 

A general problem of the change point model is that the theorem above cannot be applied due to the 
structural breaks. Even though the proposed algorithm might be a powerful tool to solve the problem, we have 
to use it carefully. Any inference after estimating the model should be backed up by some Monte-Carlo studies. 

5 . 4 . Threshold AR-ARCH model 

Threshold AR-ARCH models are popular when the mean or variance reversion properties change dependent 
on the past of the process. Threshold AR models are popular as they are simple but powerful examples for 
regime switching models. Threshold ARCH processes have many applications in finance, because they are 
suitable to capture the so called leverage effect. 

The general model is given by 

j&Vk&Iij I 


with thresholds and 


— “i.o + X/ X/ X/ > bk,i}\ej^t-k\^' + Si^t 

jevkeiij i 

with thresholds b^^i. The option of one threshold at bij. = 0 in the conditional variance model is very popular. 
This leads to the well known TARCH model, introduced by Rabemananjara and Zakoian (1993). Ziel et al. 
(2015) applied the proposed algorithm to a similar multivariate AR-TARCH type model to electricity market 
data. Here, we can use the algorithm proposed above, because all covariate processes and Yt can be weakly 
stationary. The mentioned zero-threshold option is often suitable in practice as it only doubles the volatility 
parameter space. 

5.5. AR-ARCH model with quadratic interactions 

Interaction models are very popular in classical regression settings, especially in medicine. This type of 
model was e.g. analysed by Choi et al. (2010) or Bien et al. (2013), but not in a time series context. In general 
we can apply the theorem for these models as well, as the interactions are in general weakly stationary processes, 
if they have still a finite second moment. The full quadratic interaction model is given by 

+ EE EEE E 4^i.j.k^l.mYj^t — kTl^t—m T 

j&Tikeli.j jGT> leT) kelij mGlij 

A problem that arises is the size of the parameter space which is Pn{Pn + l)/2, where the standard AR-ARCH 
model has Pn parameters. 

5.6. ARMA-CARCH model 

The last extension considers a very popular class of models. We know that every ARMA(p, q) model can 
be rewritten as an AR(oo). Similarly a univariate GARCH(p, q) can be expressed as an ARCH(oo). Hence, it 
is clear that every ARMA-GARCH model can be written as an AR(oo)-ARCH(oo). This AR(oo)-ARGH(oo) 
can be well approximated by an AR(^-ARCH(5) for large p and q. However, this gives an approximation and 
will likely include more parameters than the original ARMA-GARGH model. 

Recently, Chen and Chan (2011) proposed a method of how to estimate ARMA processes in a lasso frame¬ 
work, using this kind of approximation. The idea is simple: Given the ARMA model 

Yi,t = 4'i,0 + 4‘i,j,kYj,t-k + dij^k£j,t-k + £i,t 

jG'DfeG/i.j keKi^j 
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we consider first an AR(^-model with large p that can approximate the true ARMA model sufficiently well. 
The residuals of this fitted model are used for constructing the regressor matrix that contains the lagged 
autoregressive part and moving average part. We repeat the lasso estimation with this regressor matrix. So this 
procedure leads automatically to a two step approach. Clearly, we can iterate this more often to receive better 
stability, similarly to the algorithm we presented. Chen and Chan (2011) showed that under certain conditions 
this estimation principle based on the adaptive lasso can lead to consistent estimates. 

The same principle can be applied to the GARCH model as well. So we hrst estimate a high dimensional 
ARCH model and take the estimated conditional variances for constructing the response matrix required for the 
GARCH model. This method opens a lot of possibilities for applications in financial frameworks. In multivariate 
settings, we have to specify a special GARCH model. In fact we can use every GARCH model that we can 
express in regression form, so even the BEKK-GARCH is possible. 

6. Simulation study 

In this section we perform Monte-Carlo simulations to learn about the finite sample properties of the model 
algorithm. Of course the results of the simulation will very much depend on the true model. For illustration 
purposes we restrict ourselves to a univariate settings where both and are increasing with a rate of ^/n. 
For all simulations we consider a one-dimensional AR-ARCH-type process 

Yt= Y. c^kYt-k + et ( 10 ) 

fce/i,i 


where e* 


(TtZt with Zt ~ N{0,1 ) and 


o't — o:o + ai\£t-i\ + a2\£t-2\ 

with ao = O.OI and ai = a 2 = 0.49. The true subset Ii^i of relevant lags of model (10) is given by Ii^i = {n^|n S 
N} = {1,4,9,16, 25,...}. For parameters (j)k with k G Ii^i we define (pk = 0.95((/>“^ — l)(j)Yk (j) = 0.85. As 

— 1) J2keii 1 1) SfeeN = 1 we have X)fcG/i i “ 0.95. So the considered process has a clear 

autoregressive structure and is stationary. In Figure 1 the considered coefficient structure and some simulated 
sample paths are visualised. In the sample paths we observe the clear conditional heteroscedasticity. For the 




(a) First considered coefficients with corresponding lag. 

Figure 1: Considered parameters in la and simulated sample paths of 3 time series in lb of model (10). 


estimation the proposed superset Xij will be important as well. We consider the set = {1,2,..., [5^/n\}, 
so we have that ~ ^/n. 

Subsequently we want evaluate the estimation procedure on the full tuning parameter path. Therefore we 
estimate (10) for all A values on a given exponential grid A = {2^\g € G} where G is a equidistant grid from 
—4 to —18 of length 100. Additionally, we want to illustrate the impact of different information criteria. The 
information criteria that we consider are the Akaike information criterion (AIC), the Hannan-Quinn criterion 
(HQC) and the Bayesian information criterion (BIC). These are all special cases of the generalised information 
criterion (see e.g. Kim et al. (2012)) that is given by GIC(k„) = log{df )+KnK/n, where K represents the number 
of parameters in the model. We get the AIC, HQC and BIC by choosing either = 2, or = 21og(log(n)) 
or K„ = log(n), respectively. The volatility model is estimated by the methods explained in the section above. 
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The model order is assumed to be known. In all adaptive lasso estimation procedures we choose only the lasso 
itself, so T = 0. We simulated for n S {300,600,1200} with a Monte Carlo sample size N = 1000. 

After simulating the process, we estimate by the proposed iteratively reweighted lasso algorithm. The first 
simulation result is given in Figure 2. There we see the proportions of both the irrelevant and relevant included 



Figure 2: Proportion of irrelevant included parameters (black to red) and relevant included parameters (black 
to green) for n = 600 and A G A. 


parameters of all estimated parameters for the homoscedastic case {k = 1) and for the heteroscedastic with 
one additional replication (fc = 2) given a situation with n = 600 observations and the exponential grid tuning 
parameter grid A. Obviously, we observe that for both models the probability to include a parameter increases 
with decreasing A. We see that parameters (j)k with k G Ii.i and small k are easier to detect than those with 
larger k. This is clear as (j)k with k G Ii.i is decreasing in k. Further, we can observe that for both cases {k = 1 
and k = 2) the algorithm seems to distinguish well between relevant parameters and irrelevant parameters. In 
this situation a reasonable choice of the tuning parameter could be log(A) = —6. There we see that proportion of 
relevant parameters (green colored) included is clearly closer to 100% than the proportion of irrelevant included 
parameters (dark red to black). It seems that the heteroscedastic algorithm can distinguish better than its 
homoscedastic counterpart. 

To emphasis this fact we created a new plot where we visualise the computed mean proportion of all irrelevant 
included parameters against the mean proportion of all relevant included parameters. The mentioned plot is 
given in Figure 3. We additionally added the corresponding values for the considered information criteria. To 
understand the impact of the sample size and the number of iterations we plot the cases for n = 600 and 
n = 1200 and the first three iterations of the algorithm. The bottom left corner corresponds to very large A 
values where no parameter at all is included in the model. The top right corner covers the ordinary least square 
estimate with A = 0. 

Roughly speaking we are aiming for estimators that are as close as possible to the upper left corner. It is 
particularly important to mention that for increasing n we should get close to the upper left corner. This seems 
to be satisfied for the relevant tuning parameter path. We see that with the heteroscedastic cases with k = 2 
and fc = 3 have better selection properties than the homoscedastic case. The improvement from the case k = 2 
to A: = 3 is very small, but it is still there. The same holds for the considered information criteria. Note that 
even though it is well known that the AIC is inconsistent in parameter selection in a finite sample setting it 
seems to perform quite well. 

Nevertheless, it is not clear how the algorithm performs in an out-of-sample forecasting study. Therefore, 
we conduct another simulation study where we focus on the out-of-sample forecasting error. We compute the 
1-step ahead mean absolute forecast error (MAE) which is defined by MAE = \^n+i — Yn+i\ where 

Yn+i denotes the forecast of Yn+i- Additionally, we calculate the forecasting error for the corresponding oracle 
model. Eor the oracle we assume that the underlying lag structure of the autoregressive model is known. 

The simulation results for n = 300 and n = 600 are given in Figure 4. We see that the homoscedastic 
algorithm performs significantly worse than the heteroscedastic one with fc = 1, except for large A values in 
the n = 600 situation. Interestingly for n = 300 and k = \ the MAE hardly goes below the value of the case 
with very large A where = 0 for all A; G Ii^i. In contrast for A: > 1 an improvement in the forecasting 
performance to the case with very large A where (()/j = 0 is possible. The same fact can be observed, within the 
oracle procedures, but the improvement is not that obvious. From an applications perspective, this is extremely 
significant. It indicates that we can benefit more from taking the heteroscedasticity into account in settings 
















6 Simulation study 


12 



Figure 3: Mean proportion of all irrelevant included parameters against mean proportion of all relevant included 
parameters for the first three iterations and n S {600,1200} on the full A grid and for the considered information 
criteria. 




(a) n = 300 


(b) n = 600 


Figure 4: MAE for n = 300 (4a) and n = 600 (4b) of the iteratively reweighted lasso (IRL) method for several 
iterations k G {1,2, 3,4}, such as their oracle estimators for the AR-ARCH model. 


with unknown model structure than in a setting where the underlying structure is known. However, we usually 
do not know the true underlying model as the oracle does, especially in high-dimensional settings. It shows that 
the proposed estimation algorithm can lead to crucial improvements in a high-dimensional setting. This is also 
observed by Ziel et al. (2015) in applications of the proposed estimation algorithm to electricity market data. 

As a robustness check we replicate the simulation study with a different volatility model. We assume a 
TARCH process for the residuals. TARCH models are popular in financial applications as they are able to 
capture leverage effects. The considered TARCH process for the simulation study is parameterised through 

Ct = ao + ai|£t_i| -|- l{£t_i < 0}|£(_i| -I- a2\St-2\ + <^2 ^{^t-2 < 0}|£t_2| 

where the leverage effect is modelled by the two parameters and which give an additional impact on 
negative past residuals to the volatility. The selected parameter setting is ai = a 2 = 0.245 and = 0.49. 

We compute the 1-step ahead mean absolute forecast error (MAE) for n = 300 and n = 600. The simulation 
results with the corresponding oracles are given in Eigure 5. There we observe similar behaviour as for the AR- 
ARCH model in Eigure 4. As there is a clear improvement in the MAE it shows that the iteratively reweighted 
lasso algorithm can work well for data with asymmetric volatility. 





































8 Summary and Conclusion 


13 




(a) n = 300 


(b) n = 600 


Figure 5: MAE for n = 300 (5a) and n = 600 (5b) of the iteratively reweighted lasso (IRL) method for several 
iterations k G {1,2, 3,4}, such as their oracle estimators for the AR-TARCH model. 


7 . Applications to electricity market data and metal prices returns 

In this section we briefly show two applications of the proposed model to real data. For both applications a 
two-dimensional AR-ARCH model to the process (Yt)tei. = (^i,tj ^ 2 ,t)tGZ is considered. 

In the first application we use the hourly day-ahead electricity spot price for Germany/Austria at the 
European Power exchange (EPEX) as one process (Yi_()tgz and the hourly electricity load of Germany as 
(l 2 ,t)i 6 Z- The considered time range is from 28.09.2010 to 17.04.2014. Eor the second example we take the 
hourly intra-day returns of gold and silver prices in U.S. dollar (from London Bullion Market Association), 
denoted as XAU/USD and XAG/USD. Here represents the gold and {Y 2 ^t)tez, the silver price returns. 

The data covers 12 years of observations from 01.01.2002 to 31.12.2013. 

Note that electricity prices are known to have a strong correlation structure. In contrast, we expect either 
no or a very weak autoregressive dependency structure for the commodity returns. 

For both applications we suppose that Yt follows an AR-ARCH model as given in (4) and (5) . As the elec¬ 
tricity data has usually a long memory we propose for the autoregressive parameters the lags lij = (1,..., 700} 
for i,j G {1, 2} and similarly for the ARCH part parameters we take Jij = (1,..., 700}. This covers a mem¬ 
ory of more than 4 weeks. For the metal prices we take lij = {1,...,200} for the conditional mean and 
Ji^j = {1,..., 200} for the volatility part. The index sets are sufficiently large to capture possible weekly de¬ 
pendencies. We consider the conservative BIC as information criterion and for the adaption parameter t we 
take the lasso case with r = 0. Then we apply the iteratively reweighted algorithm and stop after i?max = 3 
iterations. Hence we solve d7?max = 2x3 lasso problems in for each application. 

The estimated /3j „ for i G (1, 2} and both applications are given in Figure 6. Here we see that in general 
most of the parameters are not included in the model. For the electricity price model there are 133 parameter 
included and for the load model 416. This matches a proportion of included parameters of 9.5% and 29.7%. 
We see that the complex autocorrelation structure that is driven by daily and weekly seasonal effects is well 
captured. 

For the metal prices we observe a different situation. HHere the gold price returns have no significant 
parameter at all. However, the silver time series exhibits a weak dependency structure. Most distinct is the first 
lag pattern with a positive coefficient for the gold returns and a negative one for the silver returns. Furthermore, 
we have two small silver coefficient clusters, a positive one around a lag of 16 hours and a negative one around 
a lag of 24. 

8. Summary and Conclusion 

An iterative algorithm to solve adaptive lasso time series problems with conditionally heteroscedastic resid¬ 
uals is described. We showed the sign consistency and asymptotic normality in a rather general time series 
setting. The asymptotic theory shows that a significant estimation improvement is possible if the conditional 
heteroscedasticity is considered. We discussed the application to AR-ARCH type models and showed applica¬ 
tions to intra-day electricity market and commodity data. 

The simulation studies underline the asymptotic results. Additionally, we showed that considering the 
heteroscedasticity in high-dimensional settings with unknown parameter specification is more important than 
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(c) Estimated coefficients for the gold price returns 


(d) Estimated coefficients for the silver price returns 


Figure 6: Estimated parameters /3j „ for the electricity market model in 6a and 6b and for the metal prices in 
6 c and 6d. 


in cases where the true underlying model is known, as it can substantially improve forecasting performance. 
This observation will likely have a strong impact on high-dimensional time series modelling, as almost every 
time series exhibits conditional heteroscedasticity, especially in economics and finance. 

The asymptotic theory shows that only two iterations are required for receiving optimal asymptotic be¬ 
haviour. Thus, the algorithm is suitable for applications, as the computational effort is only doubled in com¬ 
parison to standard homoscedastic situations. 

For future research it might be important to analyse the mentioned model extensions more carefully. Another 
very important issue is to identify the optimal penalty parameter A„ in high-dimensional time series settings. A 
different direction of further research might concern the robustness of the algorithm. The performed simulation 
study carried out that the algorithm works well in a finite sample setting. However, in a heavy tailed situation 
, it might be worth considering the LAD-lasso (see e.g. Wang et al. (2007a)), which minimises the sum of the 
absolute residuals, instead of their squares (as in lasso type algorithms). Another direction that seems to be 
a promising extension concerns the £q penalty itself. An extension to elastic net estimators, which combine ii 
and £2 penalties, could also improve estimation power. Recently Gefang (2014) applied the elastic net method 
successfully to homoscedastic multivariate AR processes. 

9. Appendix 

Proof: Theorem 1. We show the sign consistency first and then the asymptotic normality. As mentioned, the 
proof extents mainly methods from Wagener and Dette (2013). Denote e„j- the j’th unit vector in M'J", a =s b 
holds if sign(a) = sign(6) and |j • 11.,/,^ Orlicz norm with ipdix) = exp(a;‘^) — 1. In proof we will introduce at some 
points several constants Ck that are positive. 

Let fc > 1 and assume that the theorem holds for k — 1. Following the Karush-Kuhn-Tucker conditions we 
have that 

(W„ - X„/3)'(wL"-'')'(F„ - a:„/3) + Xnv'J^l 

is minimised by /3 = (/3(1)', O')' G if and only if 

A,(l)'(Tylf-il)2(F„ - X„/3) = sign(/3,) if /?, ^ 0 and 
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- X„/3)| < if 13, = 0 

holds. Thus, we have the estimator /3[f^ = (1)', O')' e K^" where 

/3L'1(1) = /3°(1) + ^(r^’(l))-'X„(l)(wL'=-'')M - ^(fL"'(l))-ia°(l) (11) 

where s°(l) = (ui,..., Wg„)'sign(/3° (1)). 

Now we define the expressions 

Vi,j = e;,,-(^’(l))-iX„(l)'(Ty[f-il)2£0 

V2,j = e(,j(r||''(i))-is°(i) 

rys., = A,(l)'(M^lf'’)"an - n-iX„(l)(f L"’ (l))-'^n(ir 

774., = A„(2r7)-iA,(l)'(Wlf-'l)2X„(l)(f^’(l))-iaO(l). 

As in Wagener and Dette (2013) we can use the argument of Huang et al. (2008) that the KKT conditions 
are satisfied if 


1^3.,- 


V4,j \ < 



( 12 ) 


holds for all j > Qn- 

Hence we receive with (11) and (12) that 


P 


(/3lf' /3°) < PiAi) + P{A 2 ) + P(A 3 ) + P(A 4 ), with 


Al = I ^ bi.il > ^bil for some j < g„| , A 2 = |^|772,, | > |/3°| for some j < qn 
A 3 = jba.il > ^fi for some j > g„| and A 4 = ||? 74 ,i| > for 


some j > qn 


So we only need to show that P{Aj) —0 as n —i 00 . 

Regarding P(Ai) we have with definition of 6 „ (see (d)) that 




< P (- max |? 7 °’°°| > +P (- ™ax - 77 ? | > +P (- ™ax 177 ? ■ - 77 °’°°! > ^ 
Vii i<i<9" 4/ \ni<j<q„ 8j \ni<3<q„ o.j g 


(13) 


where 77 ?,= e'„,^.(r° (l))-iX„(l)'(W^« ) 2 s 0 and 77 ?;“ = e'„,^.(r° (l))-iX„(l)'(W^«) 
For estimating the first term in (13) we observe that 

^e'n,,(T nW)-^ X nil)'{W^nY 
V n 


■ 012^0 


< 


< 


(rji))-' 


-^Xn{iy 


l|W^°ll2 


(r°(i))-i J|r°(i)|L^ < (r°(i)) 




for sufficiently large n with ||W °||2 < CTmin by (g). Furthermore by assumption (f) we know that 

||r°(l)||| =Op(l) and ||(r°(l))-i|| 2 =Op(l). 


11 "P^ M 11 ^ (T^ 

n Vo ^min 

2 " 


(14) 


Thus we get that 


(||e;,,-(r°(l))-iX„(l)'(T^°)i < A-^^yA^'^min) 


^ 1 


for n —i 00 . With Lemma 1 (i) of Huang et al. (2006) and tail assumption (k) we can deduce that 


0,00 


< 


ipd 


-^ey^{Tlil))-^Xn{iy{W°nrel^n < Cl log(n)if‘'=A 


(15) 




































9 Appendix 


16 


for sufficiently large n, as ||A|j 2 < c||A||^^ for some c > 0 . 

Thus, we can conclude with Markov inequality, Lemma 2.2.2 of Van der Vaart and Wellner (1996) and (15) 
that 


1 


maxi<j<q^ \ri^ 


0 ,oo I 




bnU 


I maxi<^-<,„ IVi’flUi 


4|| maxi<j<,„ \r] 


i;ri 


-1 


< 4’d 

< V'd 

< V'd 


0 ,oo I 


4||maxi<j<,^ Iffiy 


.j 

br,.n 


^4c2'ipd ^(g„)maxi<j<,„ \\Vi’]°\\'>Pd 


bnVn 


4 c 2 log(l + Qn) ^ C 2 log(n)^ffi-l} 


(16) 


as V'd ^(a^) = log(l + x)^- Hence by assumption (j) we have P maxi<j<q^ \Vi’^\ ^ “t 0. 

If assumption (k) is not satisfied we can not use equation (15) to derive that P maxi<j< 5 ^ \Vi'^\ ^ “t 
0. But we can conclude with Chebyshev’s inequality and (33) shown below that 


P 


(I 

— max 
\n i<3<q. 




16 

rPbl i<j< 9 „ 


0 ,oo |2 
3 


) = Of 



Thus even without assumption (k) it holds with assumption (j) that P ma.xi<j<q^ ^ V) How¬ 

ever, note that either (k) or (k’) is required for estimating the probability of A 3 in a similar situation. 

For the second term in (13) we proceed as in Wagener and Dette (2013). We get 




(^(r°(i))-ix„(i)'((TV °)2 


(Ty[(=-il)2) + (r°(i) _ r^’(l))X„(l)'(TV[f-il) 2 ^ £° 


<||((w°)2 - (ivif'¥)^n(i)(r°(i))-'e,-„||2|l£°||2 + |l(iv[f-ii)2x„(i)||2||£°(r°(i) 

<||(tv°) 2 - (Tvlf-'’)" Il 2 ||nr°(i)|||||(r°(i))-I|| 2 |l £°||2 

+ ||(W°)2||2||nr°(l)|||||(f°(l))-i - (fL"’(l))-i||2||eO||2. 


rL"(i))|| 


2 


(17) 


All these appearing single norms we will estimate now. 

For the estimation of ||(W°)^ — (W[f “^^)^||2 we get directly with assumption (g), (h) and (i) that 

\\{Wy-{W'^t^^fh = Op(^^y (18) 

For estimating ||£° II 2 the triangle inequality yields ||e° II 2 < ||e° - e^.nlb + Ik^.nlb- We have 


k°-e^.J |2 = 




0 


(19) 


l=Pr.-|-l ^ 

11^^,rail = Op{y/n) by law of large numbers. Thus we get 

Ikralb = Op{y/n). (20) 

Further we have ||W[f “^^||2 = Op(l) by assumption (g). 

Next, we have as in Wagener and Dette (2013) that with assumption (f), and equations (14) and (18) that 

||f°(i) - f^’(i)||2 < |ir„(i)||2||(TV°)2 - = Op . 


This leads to 


(r°)-i-(r^’)-i| 


. = Op[^ 


( 21 ) 
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by using the the triangle inequality 

\\A-^-{A + B)-^h<\\A-^-{A + B)-^+A-^BA-^h + U-^BA-^\\2<Op{\\B\\2) + \\A\\l\\Bh = Op 

for two matrices A = r„(l) with i3 = r„(l) — r„ (1) and the Taylor series expansion of {A + B)~^ 
A-\ 



around 


Using all the estimated norms ((14), (18), (20) and (21)) we receive for (17) that 

-(W"L"-'’fll 2 ||nr°(i)|!|||(r°(i))-i|| 2 ||£^ 

+ ||(W°)2||2||nrO(l)|||||(f°(l))-l - (fL'’(l))-l||2||£°||2 

<Op{^)Op{^)Op{\)Op{^) + Op(l)Op(v^)Op(-^)Op(v^) = Op{KVn). 

vn Vn 


Thus we get, 

- max \r]ij - Vi j\ = -Op{hnVn) =Op (. 
n i<j<q„ n \VnJ 

This yields with assumption (j) that P maxi<_,<q^ \r]ij — rj^ j \ > —>■ 0 as n —)■ oo. For the third term in 

(13) we get with (g), (14) and (19) that 


^\vIj - ViTl < ^|e'„,,(r„(l))-ix„(l)'«)2(£° - e°,^)| 


<ll(rji))-^||2 


\/n 


^°ll 2 lk^-<ool |2 


<||(r„(i))-i||2||r! 


nlll l|W^nll 2 ll-„ 


112^0 


as n —)■ OO. So we have that —t oo as n —)■ 0. This implies T’(Ai) — )• 0. 

Now we consider P{A 2 ) < P(^ maxi<j<,„ |ry 2 j| > &«)• We have \r/ 2 ,j\ < ||(r„ )“^|| 2 |ls° (1)||2 for each 

_Q _ 

j G {1, ■ ■ ■ , 9 ra}. By Weyl’s perturbation theorem for the matrices (r„(l))“^ and (r„ (1))”^ we have for each 
ordered pair of eigenvalues that 


A, (rji))-^ -A, (r„ (1)) 




< ||(r°(i))-i-(r^’(i))-i| 


for each j G {1,..., g„}. As ||(r„(l))“^ — (r„ ( 1))“^||2 —t 0 in probability, we get with (f) that 

ll(fi'')-'|| 2 < Arfo„ + C 3 ( 22 ) 

with probability arbitrarily close to 1 for sufficiently large n. 

Furthermore, with assumption (d) we have 


||s°(l )||2 < max lAnitjl"^ < 

Y vbn 


Hence we have with assumption (j) that < P (^f maxi<j<q^ \V 2 ,j \ > ^n) < P 
For A 3 we receive similarly as for Ai that 

P(A,) <P + P |,J, - > Y^) 

\ ^ \ / \ 


(23) 


0 as n —>■ 00 . 


+ P [ max |? 73 ,j - + pf max 

\g 7 z<J<Pn Id / \qn<J<Pn. 


(24) 


where 

= A,(1)'(TU«)2(/„ - p-ix„(i)(r°)-ix„(i)'(FF°)")£°, 
7?°;“ = X,{inwy{I^ - n-iX„(l)(f°)-iX„(l)'(MA°)^)£°, 
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As in Wagener and Dette (2013) we consider 773 ’°° = with 

= X,{l)\Wlf{Ir. - n-iX„(l)(r°)-iX„(l)'(TyO)2). 

Then we have for sufficiently large n that 

iK.lb < ||A,(l)'||2||(Ty°)^||2(l + ||n-iX„(l)(r°)-iX„(l)'|i2||(W°)2||2) 

= Op(V^)Op(l)(l + Op{l)Opil)) = Op(V^). (25) 

Now we receive we receive as in equation (15) with Huang et al. (2006) Lemma 1 (i) and assumption (k) that 




0,oo 


< 


V’d 


-H° £° 

~ 71,j^(x>,n 


< ce\og{n)^^‘^ 


(26) 


V’d 


Thus we get with Markov inequality, Lemma 2.2.2 of Van der Vaart and Wellner (1996), (26) and assumption 
(j) similarly to (16) that 


P 


max 

\<]n<j<Pn 




<1pd 

<i^d 




\8cei^d ^(9") maxg„<j<p^ hsjIUd / 

f 

Vc7v^log(l +p„ - g„)a log(n)i{‘^=i> 


(27) 


as n —)■ 00 . 

If instead of (k) the alternative assumption (k’) holds we can not use equation (26) to derive that it holds 
P ^maxq^<^<p^ 1 ^ 3 ’^ I > 0 . But can get with Chebyshev’s inequality and (33) shown below that 


max 177 ! 


qn<j<p 


A’ri ^ 




< 


64 

Kxl 


'^n n 


j = qn + l 


Thus with (k’) it holds P |^max,j^<j<p^ \Vi’j°\ ^ 0- 

For estimating the second term in (24) we note with (19) and (25) that 


0 O.oo 


< 


.(e° - e° ) 

n,i\ n 00 ,n/ 


< 


Cs 




l 2 ||£°-£° I 


= —/=C)p{^/n). 

vn 


Hence we have with assumption (j) that 


P 


max 

V'?n<i<p. 




16 J 


< P 



-X 0 


(28) 


as n —)■ 00 . Now we estimate the third term in (24). As in (17) we get the estimate 


\vh - V 3 ,j\< x,iiy - n-ix„(i)(rj-ix„(i)'(iv):)2) 

- (TV|f-'')^(4 - n-iX„(l)(fi'’)-iX„(l)'(TV[f-'l) 2 ^ £° 

<||A,(l)|| 2 ||(TV°) 2 -(IvIf-'') 2 || 2 ll £°||2 


(29) 


+ 11 ^.( 1)112 


x„(i)(r°)-ix„(i)'(TV °)2 - x„(i)(ri''Vix„(i)'(ivL '=-'')2 


lk °||2 


using estimates derived for Ai. For the lengthy norm in (29) we get 


(^x„(i)(r°)-ix„(i)'(iv°)2-x„(i)(r^yix„(i)'(iv|f-ii)2^ ^ 

<||(tv°) 2 - (TvL"-'’)"ll2||r„(i)|i2ll(rL"Vi2||(T^L"-'')"l!2 
+ ||(f°)-i - (fr'Vi2l|r„(i)||2||2ll(T^™)"ll2ll2ll(^^L'=-'')"^ 

+ Wiwlf - (TvL'-'i)2||2||r„(i)||2||(fr)-i2||(M^™)"||2 

=0p( + )0,(1)0,(1)0p(1)=0p( + ) 
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by equation (18), (22), (14), (21), 
for (29) with (20) that 


< CTinin and ||W[f ^^||2 = Op{l) by assumption (g). Thus we receive 


K, - = Op{^)Op Op(v^) + Op{^)Op Op{^) = Op{KVn) (30) 

Hence we have with assumption (j) and (30) that 




< P 




> Cio —t 0 


(31) 


as n —>■ oo. Thus we get for (24) with the estimates (27), (28), (31) and assumption (e) that PiA^) —>• 0. 
For missing event A 4 the situation is similar. We have that 

^^(^ 4 ) l^4j| > ^ (g ^ ■ 

As it holds with (23) that 


I-., I < ^ 

^77 

< — 
- 2 


^.■(i)'(w^L'-'')"^n(i)(ri''(i))-' 


N^(i)lh 


^^n(l)' 




(ri"'(i))-' 


K{i)\ 


= X^Op{l)Op{l)Op{l)Op 


\/^n 




^ny/On 

'\/^n 


we get with assumption (e) and (j) that ^^(^ 4 ) —)■ 0 as n —)■ 00 . Hence, /3jf^ is sign consistent. 

For the asymptotic normality we use similar concepts as in Wagener and Dette (2013). So given sign 
consistency of we have from equation (11) that 


n In 


(32) 


If we subtract /3° (1) and multiply the result by we directly get 


/n 


e;(/3W(i)-/3°(i)) = 


1 


(fc)‘ 


Arj 


^;(rr)-'x„(i)'(M^if-iy£o- , yy,, e;(ri"’)-^a°(i). 


Sn{k) " " "■ i/ns„(fc) 

For the second term we get with (22), (14) and H^nlb = 1 that 


2y/nsn(k) 


Xn 


2y/nsn(k) 




< 


Xn 


< 


2y/nsn{k) 
Xn^ Qnb 


2Sn{k)\/nbn 


ll&lljIKrllV'lblkSii)^ 


(Arh„+ ■=.) = Op 


With assumption (j) this converges to zero. 

For estimating the first term we use the decomposition 


(fi"’)-iX„(l)'(Wy'y = Hi + 52 + 53, where 5i = (f ° )-iX„(l)'(FF° )2, 

52 = ((rL'”'’)-' - (r°)-i)X„(l)'(M^°)2, and Ha = - Wlf). 

Now we decompose - £)(„_„). For the first term we have 

= Er=i with at = Cn(r°)~^X„,4l). So we can calculate EXlLi = 0 and 

£(^^1 E(at)^E(Z()^ = 1. It holds with assumption (j) that 

max |at| < ^ J |gn|| 2 ||(r°)~y 2 max ||cr-iX„_t(l )||2 < ^ max ||X„,t(l )||2 = O C^^'^ ) 0 

l<t<n y/nSnyk) l<i<n yjn l<t<n yjn 

for n —)■ oo. So the Lindeberg condition is satisfied and we get with the central limit theorem that 

-AnBie^ 


1 


lik) 


.0 

'oo,n 


iV( 0 ,l) 


(33) 
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in distribution as n —)■ oo. Moreover we obtain 


y/nsn{k) 




y/nSn{k) 

Cl2 

y/nSn{k) 


||(rj-i||2||x„(iri|2||(M^“)^||2||£°-£: 

v^ll^n - ^ 0 


=-0 |L 

'oo,n M ^ 


(34) 


as ||£° - eSo,nll 2 ^ 0 as n ^ oo. 

Regarding B 2 we similarly to Wagener and Dette (2013) that 


fnsn{k) 


i'nB2el 


< 


< 


r ..J ICn||2||((ri" ")-' - (r°)-i)x„(i)'(TyO)2£0||2 

^nSnK,k) 


O'minV ^ 

Using triangle inequality we get 

||X„(1)'(W°)2£°II 2 < \\Xr,{l)'{Wlf{el - £^,J||2 + \\Xr.{iy 

For the first term we have as above 

\\Xr.{l)'{Wlf{el - £^,J ||2 < \\nT°J-i\\W°M\el - = OpiV^)Op{l)Op{l) = Op{V^). 

For the second term we get with Markov’s inequality 


(35) 


P 


( _}_\W (■\xr0\2^0 ||2 


\qnn 


ix^iinw'iye^. 


>c)<—E e(e^m-) <4- 




for c > 0. This gives ||X„(1)'(W°)2e;^ ^||2 = 0{y^). With IKril" ^') ^ - (F^) = C)(^) and the 

previous estimates it follows for (35) that 




/nsn{k) 


enB2e° 


= 0[^]0[^]0 (v/^) = O 


which converges to 0 with assumption (j). 

For the last term that corresponds to we have with (j) that 


yk) 




< 


\Acv 


(ApLn + C4)||X„(1)'((TU°)2 - iwy-^^r)slh 


^n -\/Qn \ 


(36) 


Again, the second norm can be estimated by 

\\XMyiiW°S - {wy-^Y)elh 

< ||X„(1)'((W°)2 - - £°^,J||2 + ||X„(1)'((TU°)2 - {Wy-^^f)el 

using the triangle inequality. The first term can be estimated by 

|1X„(1)'((FF°)2 - (WE'l)2)(e° - e^,J||2 < ||X„(l)'|i2||(T^O)^ - 

< Op{Vn)Op Op{l) = Op{K) 

For the second term we have again with assumptions (h), (i) and and Markov’s inequality 




q-n / n 

P (||X„( 1 )'((W °)2 - > c) < E®: E 


Z=1 \i = l 


- (ay-y 


-£t 


^ i^l \t^l J 

where up for 1 < t < n are the diagonal elements of and c > 0. Hence with assumption (j)we 


receive 


_i_ tf R c-O I / ^ -/^hn 

/7\4'^3^nl — ^8 /— ^ ^ 

Sn\tx) yU 


(37) 
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as n —>■ oo. With the three estimates involving Bi, B 2 and B 3 we receive for equation (32) together with 
equations (33), (34), (36), (37) and Slutky’s theorem that —/3°(1)) —)■ -/V(0,1). 

At the beginning that the theorem is satisfied for fc > 1. So the proof of the inital step with fc = 1 is 
missing. However, the proof is similar to the sign consistency and asymptotic normality proof with fc > 1, as 
Wagener and Dette (2013) explained it for the unconstrained weighted adaptive lasso. Note that the proof itself 
is less complex than the case fc > 1, but involves the eigenvalue assumptions to the unsealed Gramian r° (i.e. 

Ao.min < Amin(r°)) that Were not used in the previous part, instead of the assumption to the scaled version r„. 

□ 
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